Import Python modules:

  • Shapefile to read shp files
  • pandas to read mobility data
  • numpy for array manipulation
  • bokeh for interactive plotting
In [1]:
import shapefile

import numpy as np
import math

import pandas as pd

import datetime
from datetime import date

from bokeh.models.widgets import DateSlider, Slider
from bokeh.models.callbacks import CustomJS
from bokeh.layouts import column, row
from bokeh.models import Glyph, Range1d, ColorBar, NumeralTickFormatter, HoverTool, LinearColorMapper, TapTool, Patch
from bokeh.plotting import figure, show, ColumnDataSource, output_file
from bokeh.palettes import RdYlGn11
from bokeh.io import output_notebook
output_notebook()
Loading BokehJS ...

Read in United Kingdom shapefile (second level subdivisions) downloaded from [GADM].(https://gadm.org/download_country_v3.html)

In [2]:
boundaryshp = "shapefiles/gadm36_GBR_shp/gadm36_GBR_2.shp"

sf = shapefile.Reader(boundaryshp)

shapes = sf.shapes()
records = sf.records()

with shapefile.Reader(boundaryshp) as shp:
    print(shp)
shapefile Reader
    183 shapes (type 'POLYGON')
    183 records (14 fields)

Records file contains the name of the corresponding region shape and additional labels.

In [3]:
records[0]
Out[3]:
Record #0: ['GBR', 'United Kingdom', 'GBR.1_1', 'England', '', 'GBR.1.1_1', 'Barnsley', '', '', 'Metropolitan Borough', 'Metropolitan Borough', '', 'GB.BX']

Convert the shapes and records into (x,y)-coordinates to be plotted as patches using bokeh.

  • region_xcoords: x-coords for each region boundary
  • region_ycoords: y-coords for each region boundary
  • region_names: names for each region boundary
  • region_codes: Area_Codes for each (sub)region boundary
In [4]:
region_xcoords = []
region_ycoords = []
region_names = []
region_codes = []

for index, shape in enumerate(shapes):
    record = records[index]

    name = record[6]
    code = record[-1]
    
    offsets = list(shape.parts[:])
    
    p_start = offsets.pop(0)
    while len(offsets):
        p_end = offsets.pop(0)
        region_xcoords.append([point[0] for point in shape.points[p_start:p_end]])
        region_ycoords.append([point[1] for point in shape.points[p_start:p_end]])      
        region_names.append(name)
        region_codes.append(code)
        p_start = p_end
        
    region_xcoords.append([point[0] for point in shape.points[p_start:]])
    region_ycoords.append([point[1] for point in shape.points[p_start:]])
    region_names.append(name)
    region_codes.append(code)

Import the Global Mobility Data, and restrict to the UK data.

In [5]:
mobility_data = "mobility_data/Global_Mobility_Report.csv"
df = pd.read_csv(mobility_data)

UKdf = df.loc[df['country_region_code'] == 'GB']
/home/ollie/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3058: DtypeWarning: Columns (3) have mixed types. Specify dtype option on import or set low_memory=False.
  interactivity=interactivity, compiler=compiler, result=result)

Rename and group regions to correspond with the naming and grouping between the 2 data sets. For example we group the shapefile regions Barnsley, Doncaster, Rotherham and Sheffield to the South Yorkshire mobility data entry.

In [6]:
#  Regrouping

region_names = list(map(lambda x: x.replace("Barnsley","South Yorkshire"), region_names))
region_names = list(map(lambda x: x.replace("Doncaster","South Yorkshire"), region_names))
region_names = list(map(lambda x: x.replace("Rotherham","South Yorkshire"), region_names))
region_names = list(map(lambda x: x.replace("Sheffield","South Yorkshire"), region_names))

region_names = list(map(lambda x: x.replace("Bradford","West Yorkshire"), region_names))
region_names = list(map(lambda x: x.replace("Calderdale","West Yorkshire"), region_names))
region_names = list(map(lambda x: x.replace("Kirklees","West Yorkshire"), region_names))
region_names = list(map(lambda x: x.replace("Leeds","West Yorkshire"), region_names))
region_names = list(map(lambda x: x.replace("Wakefield","West Yorkshire"), region_names))

region_names = list(map(lambda x: x.replace("Wirral","Merseyside"), region_names))
region_names = list(map(lambda x: x.replace("Sefton","Merseyside"), region_names))
region_names = list(map(lambda x: x.replace("Knowsley","Merseyside"), region_names))
region_names = list(map(lambda x: x.replace("Saint Helens","Merseyside"), region_names))

region_names = list(map(lambda x: x.replace("Manchester","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Salford","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Bury","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Bolton","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Oldham","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Wigan","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Trafford","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Bury","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Rochdale","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Stockport","Greater Manchester"), region_names))
region_names = list(map(lambda x: x.replace("Tameside","Greater Manchester"), region_names))  

region_names = list(map(lambda x: x.replace("Coventry","West Midlands"), region_names))
region_names = list(map(lambda x: x.replace("Solihull","West Midlands"), region_names))
region_names = list(map(lambda x: x.replace("Birmingham","West Midlands"), region_names))
region_names = list(map(lambda x: x.replace("Sandwell","West Midlands"), region_names))
region_names = list(map(lambda x: x.replace("Walsall","West Midlands"), region_names))
region_names = list(map(lambda x: x.replace("Wolverhampton","West Midlands"), region_names))
region_names = list(map(lambda x: x.replace("Dudley","West Midlands"), region_names))

region_names = list(map(lambda x: x.replace("Newcastle upon Tyne","Tyne and Wear"), region_names))
region_names = list(map(lambda x: x.replace("Sunderland","Tyne and Wear"), region_names))
region_names = list(map(lambda x: x.replace("North Tyneside","Tyne and Wear"), region_names))
region_names = list(map(lambda x: x.replace("South Tyneside","Tyne and Wear"), region_names))

# Renaming

region_names = list(map(lambda x: x.replace("Poole","Dorset"), region_names))
region_names = list(map(lambda x: x.replace("Bournemouth","Dorset"), region_names))
region_names = list(map(lambda x: x.replace("Telford and Wrekin","Shropshire"), region_names))
region_names = list(map(lambda x: x.replace("Halton","Borough of Halton"), region_names))  
region_names = list(map(lambda x: x.replace("Caerphilly","Caerphilly County Borough"), region_names))
region_names = list(map(lambda x: x.replace("Torfaen","Torfaen Principal Area"), region_names))
region_names = list(map(lambda x: x.replace("South Ayrshire","South Ayrshire Council"), region_names))
region_names = list(map(lambda x: x.replace("North Ayrshire","North Ayrshire Council"), region_names))
region_names = list(map(lambda x: x.replace("Perthshire and Kinross","Perth and Kinross"), region_names))
region_names = list(map(lambda x: x.replace("Neath Port Talbot","Neath Port Talbot Principle Area"), region_names))
region_names = list(map(lambda x: x.replace("Rhondda, Cynon, Taff","Rhondda Cynon Taff"), region_names))
region_names = list(map(lambda x: x.replace("Merthyr Tydfil","Merthyr Tydfil County Borough"), region_names))
region_names = list(map(lambda x: x.replace("Highland","Highland Council"), region_names))
region_names = list(map(lambda x: x.replace("Durham","County Durham"), region_names))
region_names = list(map(lambda x: x.replace("Eilean Siar","Na h-Eileanan an Iar"), region_names))
region_names = list(map(lambda x: x.replace("Argyll and Bute","Argyll and Bute Council"), region_names))
region_names = list(map(lambda x: x.replace("East Lothian","East Lothian Council"), region_names))
region_names = list(map(lambda x: x.replace("Angus","Angus Council"), region_names))
region_names = list(map(lambda x: x.replace("East Ayrshire","East Ayrshire Council"), region_names))
region_names = list(map(lambda x: x.replace("Dundee","Dundee City Council"), region_names))
region_names = list(map(lambda x: x.replace("Glasgow","Glasgow City"), region_names))
region_names = list(map(lambda x: x.replace("Orkney Islands","Orkney"), region_names))
region_names = list(map(lambda x: x.replace("Aberdeen","Aberdeen City"), region_names))
region_names = list(map(lambda x: x.replace("Aberdeen Cityshire","Aberdeenshire"), region_names))
region_names = list(map(lambda x: x.replace("East Renfrewshire","East Renfrewshire Council"), region_names))
region_names = list(map(lambda x: x.replace("East Dunbartonshire","East Dunbartonshire Council"), region_names))
region_names = list(map(lambda x: x.replace("West Dunbartonshire","West Dunbartonshire Council"), region_names))
region_names = list(map(lambda x: x.replace("Armagh, Banbridge and Craigavon","Armagh City, Banbridge and Craigavon"), region_names))
region_names = list(map(lambda x: x.replace("North Down and Ards","Ards and North Down"), region_names))
region_names = list(map(lambda x: x.replace("Gateshead","Tyne and Wear"), region_names))
region_names = list(map(lambda x: x.replace("Anglesey","Isle of Anglesey"), region_names))
region_names = list(map(lambda x: x.replace("Conwy","Conwy Principal Area"), region_names))
region_names = list(map(lambda x: x.replace("Wrexham","Wrexham Principal Area"), region_names))
region_names = list(map(lambda x: x.replace("Bedfordshire","Bedford"), region_names))
region_names = list(map(lambda x: x.replace("Central Bedford","Central Bedfordshire"), region_names))
region_names = list(map(lambda x: x.replace("Bridgend","Bridgend County Borough"), region_names))
region_names = list(map(lambda x: x.replace("Bristol","Bristol City"), region_names))

Let's arbitrarily choose to visualise the grocery_and_pharmacy_percent_change_from_baseline data.

In [7]:
UKdf.head()
Out[7]:
country_region_code country_region sub_region_1 sub_region_2 iso_3166_2_code census_fips_code date retail_and_recreation_percent_change_from_baseline grocery_and_pharmacy_percent_change_from_baseline parks_percent_change_from_baseline transit_stations_percent_change_from_baseline workplaces_percent_change_from_baseline residential_percent_change_from_baseline
74231 GB United Kingdom NaN NaN NaN NaN 2020-02-15 -12.0 -7.0 -35.0 -12.0 -4.0 2.0
74232 GB United Kingdom NaN NaN NaN NaN 2020-02-16 -7.0 -6.0 -28.0 -7.0 -3.0 1.0
74233 GB United Kingdom NaN NaN NaN NaN 2020-02-17 10.0 1.0 24.0 -2.0 -14.0 2.0
74234 GB United Kingdom NaN NaN NaN NaN 2020-02-18 7.0 -1.0 20.0 -3.0 -14.0 2.0
74235 GB United Kingdom NaN NaN NaN NaN 2020-02-19 6.0 -2.0 8.0 -4.0 -14.0 3.0

We load in the grocery and pharmacy data replacing missing data with the string 'null'.

In [8]:
grocery_and_pharmacy_percent_change_from_baseline = []

for region in region_names:
    region_data = UKdf.loc[df['sub_region_1'] == region, 'grocery_and_pharmacy_percent_change_from_baseline']
    if region_data.empty:
        region_data = None
    else:
        region_data = np.array(region_data)
    grocery_and_pharmacy_percent_change_from_baseline.append(region_data)
    

mobility_data = []

for entry in grocery_and_pharmacy_percent_change_from_baseline:
    if entry is None:
        mobility_data.append('null')
    elif math.isnan(entry[0]):
        mobility_data.append('null')
    else:
        mobility_data.append(entry[0])

Initialise the dates over which we plot our data.

In [9]:
start_date = date(2020, 2, 15)
last_date = date(2020, 7, 19)

date_list = [start_date + datetime.timedelta(days = i) for i in range((last_date-start_date).days)]

We will plot a UK map showing the percentage change from baseline of grocery and pharmacy travel. For more details about the data see Community Mobility Reports.

We will provide a date slider to interact with the plot to see how this varies thorughout the lockdown period. Tapping on a region will plot the percentage change from baseline of that region as a line graph.

Define ColumnDataSources for the plots.

In [36]:
source = ColumnDataSource(data={'x': region_xcoords, 'y': region_ycoords,
                                'name': region_names, 'code': region_codes,
                                'colour_data': mobility_data
                               }
                         )



local_source = ColumnDataSource(data={'x': [1581724800000 + i*1000*60*60*24 for i in range(len(date_list))], 
                                      'y': ['null' for i in range(len(date_list))]
                               }
                         )

vline_source = ColumnDataSource(data=dict(c=[1581724800000],
                                          y= [-300],
                                       length=[0],
                                       angle=[np.pi / 2]
                                       )
                             )

We use javascript rather than python callbacks so we can host the interactive map on github pages.

In [37]:
# Set high and low percentage change to normalise colorbar
low = -75
high = 75
mapper = LinearColorMapper( palette=RdYlGn11[::-1], low=low, high=high)
mapper.nan_color = 'whitesmoke' # missing entries will be coloured whitesmoke
In [38]:
# Map Plot
TOOLS = "pan,wheel_zoom,box_zoom,reset,hover"

map_plot = figure(title="UK Covid Mobility Data - Change in Grocery and Pharmacy Visits from 6 Week Baseline", 
                  tools=TOOLS, toolbar_location = "below",
                  x_axis_location=None, 
                  x_range=Range1d(-10.5, 4.5, bounds="auto"),
                  y_axis_location=None, 
                  y_range=Range1d(49.5, 61.5, bounds="auto"),
                  background_fill_color = 'paleturquoise')

map_plot.grid.visible = False

color_bar = ColorBar(color_mapper=mapper,
                     label_standoff=12, border_line_color=None, location=(-20,0))

map_plot.add_layout(color_bar, 'right')

hover = map_plot.select_one(HoverTool)
hover.point_policy = "follow_mouse"

hover.tooltips = [("Region", "@name"),
                  ("Area Code", "@code"),
                  ("Percent Change", "@colour_data%"),
                  ]
In [39]:
# Local Data Line plot

local_trend_plot = figure(title='Regional Plot', toolbar_location = None,
                          x_axis_type="datetime", height = 300, width =600,
                          x_range=Range1d(1581724800000, 1595116800000, bounds="auto"), 
                          y_range=Range1d(low, high, bounds="auto"), margin= (0,0,0,0),
                          background_fill_color = 'paleturquoise'
                         )

local_trend_plot.xgrid.visible = False
local_trend_plot.ygrid.grid_line_color = 'navy'
local_trend_plot.ygrid.grid_line_alpha = 0.5

local_trend = local_trend_plot.line('x','y', source = local_source, color='navy', line_width=3)

regions = map_plot.patches('x', 'y', source=source, 
                 fill_color={'field': 'colour_data', 'transform': mapper},
                 line_color='white', line_width=1, selection_line_color = 'navy',
                           
                           nonselection_fill_alpha=1, nonselection_line_alpha=1)

local_trend_plot.ray(x="c", y="y", length="length", angle="angle",
                 source=vline_source, color='navy', line_width=2, alpha=0.5)

plot_region_data = CustomJS(args = dict(local_source = local_source, 
                                        title = local_trend_plot.title,
                                        region_names= region_names,
                                        source = source, 
                                        mobility_data = grocery_and_pharmacy_percent_change_from_baseline), 
                                        code = """

const region_index = source.selected.indices

var y = local_source.data['y']

title.text = region_names[region_index]

if (mobility_data[region_index] != null){
    for (var i = 0; i <= mobility_data[region_index].length; i++) {
        y[i] = mobility_data[region_index][i]    
    }
    local_source.change.emit();

}


""" )

map_plot.add_tools(TapTool(callback = plot_region_data, renderers = [regions]))
In [40]:
# Date Slider

date_slider = DateSlider(title="Date", 
                         start=start_date, end=last_date, value=date(2020, 2, 16), 
                         step=1000*60*60*24,
                        width = 560,
                        margin= (0,0,0,35))

date_slider.bar_color = "paleturquoise"

date_slider_callback = CustomJS(args=dict(source=source,
                                          mobility_data = grocery_and_pharmacy_percent_change_from_baseline, 
                                          vline = vline_source), 
                                code="""
var date = cb_obj.value

var time_index = Math.round((date-1581724800000)/(1000*60*60*24))

var data = source.data

var colour_data = data['colour_data']

var x_pos = vline.data['c']

x_pos[0] = cb_obj.value

vline.change.emit()

for (var i = 0; i < colour_data.length; i++) {
            if (mobility_data[i] != null)
            colour_data[i] = mobility_data[i][time_index]
        }
        source.change.emit();

""")

date_slider.js_on_change('value', date_slider_callback)
In [41]:
show(row(map_plot,column(local_trend_plot, date_slider)))

We observe the general trend that the use of grocery and pharmacy shops increases above the baseline rate around $\sim$18 March before reducing significantly during lockdown.

There are a number of further developments to add to this tool:

  • Optimise performance
    • Decrease level of detail in plot whilst interacting
    • Use alternative data structures for more efficient updates
  • Add dropdown menu to switch between each data set
  • Overlay travel information on map
In [ ]: